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ABSTRACT 

We present results from axisymmetric, time- dependent 
magnetohydrodynamic (MHD) simulations of the collapsar model for 
gamma-ray bursts. We begin the simulations after the 1.7 M iron core of a 
25 M Q presupernova star has collapsed and study the ensuing accretion of the 
7 M Q helium envelope onto the central black hole formed by the collapsed iron 
core. We consider a spherically symmetric progenitor model, but with spherical 
symmetry broken by the introduction of a small, latitude-dependent angular 
momentum and a weak radial magnetic field. Our MHD simulations include 
a realistic equation of state, neutrino cooling, photodisintegration of helium, 
and resistive heating. Our main conclusion is that, within the collapsar model, 
MHD effects alone are able to launch, accelerate and sustain a strong polar 
outflow. We also find that the outflow is Poynting flux-dominated, and note 
that this provides favorable initial conditions for the subsequent production of a 
baryon-poor fireball. 

Subject headings: accretion, accretion disks - gamma rays: bursts - methods: 
numerical - MHD - stars: winds, outflows 
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1. 



Introduction 



The collapsar model is one of most promising scenarios to explain the huge release of 
energy in a matter of seconds, associated with gamma-ray bursts (GRBs; Woosley 1993; 
Paczyhski, 1998; MacFadyen & Woosley 1999, hereafter MW; Popham, Woosley & Fryer 
1999; MacFadyen, Woosley & Heger 2001). In this model, the collapsed iron core of a 
massive star accretes gas at a high rate (~ 1M s" 1 ) producing a large neutrino flux, a 
powerful outflow, and a GRB. Despite many years of intensive theoretical studies of these 
events, basic properties of the central engine are uncertain. In part, this is because previous 
numerical studies of the collapsar model did not explicitly include magnetic fields, although 
they are commonly accepted as a key element of accretion flows and outflows. 

In this letter we present a study of the time evolution of magnetohydrodynamic (MHD) 
flows in the collapsar model. This study is an extension of existing models of MHD accretion 
flows onto a black hole (BH; Proga & Begelman 2003, PB03 hereafter). In particular, we 
include a realistic equation of state (EOS), photodisintegration of bound nuclei and cooling 
due to neutrino emission. Our study is also an extension of MWs collapsar simulations, as 
we consider very similar neutrino physics and initial conditions but solve MHD instead of 
hydrodynamical equations. 



To calculate the structure and evolution of an accreting flow, we solve the equations of 



2. 



Method 



MHD: 
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where p is the mass density, P is the total gas pressure plus radiation pressure, v is the 
fluid velocity, e is the internal energy density, $ is the gravitational potential, B is the 
magnetic field vector, J is the current density, r\ T is an anomalous resistivity, and C is the 
cooling rate due to neutrinos. 

To compute resistivity, we follow Stone & Pringle (2001, see their equations 5 and 
Al). We perform simulations using the pseudo-Newtonian potential of the central mass 
<fr pw = GM/(r — R s ), where R s = 2GM/c 2 is the Schwarzschild radius, introduced by 
Paczyhski & Wiita (1980). We increase the mass of the BH during the calculation by the 
amount of baryonic rest mass accreted through the inner boundary. 

Our calculations are performed in spherical polar coordinates (r,0,<f>). We 
assume axial symmetry about the rotational axis of the accretion flow {9 = 0° 
and 180°). The computational domain is defined to occupy the radial range 
r { = 1.5 R s < r < r = 1000 R s , and the angular range 0° < 9 < 180°. The 
r — 9 domain is discretized on a non-uniform grid as in PB03, which yields Ar/r = 0.278 at 
the inner edge of the simulations. 

We adopt a realistic EOS, which includes contributions from an ideal gas of nuclei, 
radiation, and electrons and positrons with arbitrary degrees of relativity and degeneracy 
(Blinnikov, Dunina-Barkovskaya & Nadyozhin 1996). To compute the cooling rate, we 
follow Itoh et al. (1989, 1990), taking into account thermal neutrino emission processes 
dominated by pair annihilation, as well as the capture of pairs on free nucleons. 

Our calculations use the ZEUS-2D code described by Stone & Norman (1992a,b). We 
have extended the code to include the realistic EOS, artificial resistivity, photodisintegration 
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and neutrino cooling. We included the terms due to the resistivity and cooling in an 
operator split fashion separately from the rest of the dynamical equations. For stability, 
these terms must be integrated using the time step computed based on the resistive and 
cooling time scales, respectively. We subcycle whenever either of these two time steps 
is smaller than the time step used in the MHD equations (see, e.g., Stone, Pringle, and 
Begelman 1999). Inclusion of a non-adiabatic EOS requires iterating over temperature, T 
(see MW for details). 

By including the neutrino cooling and realistic EOS, we consider very similar 
microphysics to that used in MW. Our simulations differ from those in MW in that we 
use the ZEUS MHD code whereas MW used the PROMETHEUS hydrodynamics code 
(Fryxell, Miiller & Arnett 1989). This means that we can self-consistently calculate 
turbulent stresses generated by the magnetorotational instability (MRI) and thus include 
the outward transport of energy and angular momentum (e.g., Balbus & Hawley 1991). MW 
implemented the effects of viscosity in the disk using an alpha viscosity as prescribed by 
Shakura & Sunyaev (1973). We do not consider self-gravity and nuclear burning; however, 
as in MW our simulations track the photodisintegration of helium. Our simulations span 
the radial direction from 1.5 to 1000 Rs whereas MWs simulations span from 9.5 to 9500 
R s - Thus we can capture the innermost part of the flow near the BH but follow the 
evolution out to a smaller radius than MWs simulations. 

We adopt PB03's boundary conditions and initial conditions for the magnetic field. 
In particular, the initial magnetic field is purely radial and weak (j3 = 8nP/B 2 1 
everywhere). For the initial conditions of the fluid variables, we follow MW and adopt the 
stellar model for the helium core of a 25 M presupernova star (model S251S7B@14233 in 
Woosley & Weaver 1995). The masses of the helium and iron core derived from this star are 
7.23 and 1.70 M , respectively. Similarly to MWs simulations, our simulations start after 
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the entire iron core is assumed to have collapsed to form a BH (with R s = 4.957 x 10 5 cm), 
but before the helium envelope has collapsed. The model predicts an inner radius of the 
helium envelope, Rn e , of 2.11x10 s cm. Outside this radius, we adopt the radial velocity 
from the stellar model, set vq — and assume a non-zero I. The angular momentum 
distribution is chosen such that the ratio between the centrifugal force and the component 
of gravity perpendicular to the rotational axis is 0.02 at all angles and radii, except where 
this prescription results in I > l rnax = 10 17 cm 2 s -1 ; then we set I = l max - Inside Rn e , we set 
ve = v<t> = and assume a free-fall radial velocity. We compute the initial density inside the 
helium envelope using the mass continuity equation and the mass accretion rate from the 
stellar model at Ru c , where p = 1.16 x 10 7 g cm -3 and v r = —8.81 x 10 7 cm s -1 . 

3. Results 

With these assumptions, there is only one free parameter which defines the strength of 
the initial magnetic field. In this letter, we present results from a single model for which 
P = 10 6 at the outer boundary. 

We find that after a transient episode of infall, lasting 0.13 s, the gas with / ^ 2Rsc 
piles up outside the black hole and forms a thick torus bounded by a centrifugal barrier 
near the rotation axis. Soon after the torus forms (i.e., a couple of orbits at r = r,), the 
magnetic field is amplified by both MRI and shear. We have verified that most of the inner 
torus is unstable to MRI, and that our simulations have enough resolution to resolve, albeit 
marginally, the fastest growing MRI mode. The torus starts evolving rapidly and accretes 
onto the black hole. Another important effect of magnetic fields is that the torus produces 
a magnetized corona and an outflow. The presence of the corona and outflow is essential to 
the evolution of the inner flow at all times and the entire flow close to the rotational axis 
during the latter phase of the evolution. We find that the outflow very quickly becomes 
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sufficiently strong to overcome supersonically infalling gas (the radial Mach number in the 
polar funnel near the inner radius is ~ 5) and makes its way outward, reaching the outer 
boundary at t — 0.25 s. Due to limited computing time, our simulations were stopped at 
t = 0.28215 s, which corresponds to 6705 orbits of the flow near the inner boundary. We 
expect the accretion to continue much longer, roughly the collapse timescale of the Helium 
core (~ 10 s), as in MW. 

Figure 1 shows the time evolution of the mass accretion rate through the inner 
boundary (left panel), total magnetic energy (second left panel), neutrino luminosity (third 
left panel) and radial Poynting and kinetic flux along the polar axis at r = 190 Rs (right 
panel). Unless otherwise stated, all quantities in this paper are in cgs units. Initially, 
during a precollapse phase, M a stays nearly constant at the level of ~ 5 x 10 32 g s _1 . 
During this phase the zero-/ gas inside the initial helium envelope is accreted. Around 
t — 0.13 s, M a rises sharply as the gas from the initial helium envelope reaches the inner 
boundary. However, this gas has non-zero I and a rotational supported torus and its corona 
and outflow form, causing a drop in M a after it reaches a maximum of 2 x 10 33 g s _1 at 
t = 0.145 s. The accretion rate reaches a minimum of 6 x 10 31 g s^ 1 at t ~ 0.182 s and then 
fluctuates with a clear long-term increase. This increase is caused by the contribution from 
gas with / < 2R s c, which is directly accreted (without need to transport /) from outside 
the main body of the torus. The total mass and angular momentum accreted onto the BH 
during our simulation (0.3 s) are 0.1 M and 3 x 10 39 g cm 2 s -1 , respectively. 

The time evolution of the total magnetic energy (integrated over the entire 
computational domain) is characteristic of weakly magnetized rotating accretion flows. 
Most of the magnetic energy is due to the toroidal component of the magnetic field. We 
note a huge increase of the toroidal magnetic field coinciding with the formation and 
development of the torus. Both and B e are practically zero during the precollapse 
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phase of the evolution. But at t — 0.14 s the total energy in equals that in B r and just 
0.025 s later the energy is higher than the B r energy by a factor of 50. At the end of 
simulations the total kinetic energy from the radial, latitudinal and rotational motion are 
4 x 10 50 , 6.5 x 10 49 , and 2.3 x 10 51 erg, respectively. These gross properties indicate that 
the magnetic energy is large enough to play an important role in the flow dynamics. 

The time evolution of the neutrino luminosity, L u , shows that the neutrino emission 
stays at a relatively constant level of 3 x 10 52 erg s _1 after the torus forms. We compute 
L u under the assumption that all the gas in the model is optically thin to neutrinos, so 
that L v is volume integrated C over the entire computational domain. We note that L v 
is dominated by the neutrino emission due to pair capture on free nucleons (the so-called 
URCA cooling). 

The last panel in Fig. 1 shows the area-integrated radial fluxes of magnetic and kinetic 
energy at r = 190 Rs inside the polar outflow. The outflow is Poynting flux-dominated, 
with the Poynting flux exceeding the kinetic energy flux by up to an order of magnitude. 

Our analysis of the inner flow shows that the outflow is magnetically driven from 
the torus. Soon after the torus forms, the magnetic field very quickly deviates from its 
initial radial configuration due to MRI and shear. This leads to fast growth of the toroidal 
magnetic field as field lines wind up due to the differential rotation. As a result the toroidal 
field dominates over the poloidal field and the gradient of the former drives an outflow. 
Figure 2 shows the flow pattern of the inner part of the flow at t — 0.2735 s. The left 
and right panel show density and \B^\ maps, respectively. The maps are overlaid by the 
direction of the poloidal velocity. The polar regions of low density and high B<p coincide 
with the region of an outflow. We note also that during the latter phase of the evolution 
not all of the material in the outflow originated in the innermost part of the torus - a 
significant part of the outflow is "peeled off" the infalling gas at large radii by the magnetic 
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pressure. Figure 2 illustrates that the inner torus and its corona and outflow cannot always 
prevent the low-/ gas from reaching the BH. Even the magnetic field cannot do it if the 
density of the incoming gas is too high. Therefore, we find that the outflow of the magnetic 
energy (mostly toroidal) from the innermost part of the flow does not always correspond 
to an outflow of gas (in other words, the Poynting flux and kinetic energy flux can be in 
opposite directions). 

Figure 3 shows the radial profiles of several quantities in our run, angle-averaged over 
a small wedge near the equator (between = 86° and 94°), and time- averaged over 50 data 
files covering a period at the end of the simulations (from 0.2629 s through 0.2818 s). We 
indicate the location of the last stable circular orbit by the vertical dotted line in each 
panel. 

The profiles of each variable are not simple power-laws but are rather complex. In 
particular, the density has a prominent maximum of 4 x 10 11 g cm~ 3 at r = 5R$- The gas 
plus radiation pressure is higher than the magnetic pressure. The rotational velocity is 
nearly Keplerian inside r = 6Rs and sub-Keplerian outside this radius. 

We measure the Reynolds stress, a gas =< pv r bv^ > /P, and the Maxwell stress 
normalized to the magnetic pressure, a mag =< 2B r B < f ) /B 2 >. Note that Figure 3 shows only 
the magnitude, not the sign, of the normalized stresses. We find that except for r < 2.5Rs 
and lORs ~ r ^ 12i?s, the Maxwell stress dominates over the Reynolds stress in the inner 
flow. The last panel in Figure 3 shows that the toroidal component of the magnetic field is 
dominant for r < 50 Rs- 

We have compared the cooling time scale and the advection time scale in the flow and 
found that overall the flow is advection-dominated expect for a small region inside the torus 
where the density reaches its maximum. 



4. Discussion and Conclusions 



We have performed time-dependent two-dimensional MHD simulations of the collapsar 
model. Our simulations show that: 1) soon after the rotationally supported torus forms, 
the magnetic field very quickly starts deviating from purely radial due to MRI and shear. 
This leads to fast growth of the toroidal magnetic field as field lines wind up due to the 
torus rotation; 2) The toroidal field dominates over the poloidal field and the gradient of 
the former drives a polar outflow against supersonically accreting gas through the polar 
funnel; 3) The polar outflow is Poynting flux-dominated; 4) The polar outflow reaches the 
outer boundary of the computational domain (5 x 10 s cm) with an expansion velocity of 0.2 
c; 6) The polar outflow is in a form of a relatively narrow jet (when the jet breaks through 
the outer boundary its half opening angle is 5°); 7) Most of the energy released during the 
accretion is in neutrinos, L v = 2 x 10 52 erg s _1 . Therefore it is likely that neutrino driving 
can increase the outflow energy (e.g., Fryer & Meszaros 2003 and references therein). 

Our simulations explore a relatively conservative case where we allow for neutrino 
emission but do not allow for the emitted neutrinos to interact with the gas or annihilate. 
The only sources of nonadiabatic heating in our simulations are the artificial viscosity and 
resistivity. 

Our main conclusion is that, within the collapsar model, MHD effects are able to 
launch, accelerate and sustain a strong polar outflow. We believe that this conclusion 
will turn out to be largely independent of the initial magnetic field strength in the stellar 
core, because MRI can rapidly amplify weak fields until they are strong enough to drive a 
powerful outflow. Since our simulations are non-relativistic, and cover only the innermost 
region of the collapsing star, we cannot determine whether our outflows are sufficient to 
yield a GRB. Additional driving could also be necessary. We also find that the outflow 
is Poynting flux-dominated, and note that this provides favorable initial conditions for 
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the subsequent production of a baryon-poor fireball [e.g., Fuller, Pruet & Abazajian 
(2000); Beloborodov (2003); Vlahakis & Kdnigl (2003); Meszaros (2002)], or a magnetically 
dominated "cold fireball" [Lyutikov & Blandford (2002)]. 
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Figure Captions 

Figure 1 - The time evolution of the mass accretion rate (left panel), total magnetic energy- 
due to each of the three field components (second left panel), neutrino luminosity (third 
left panel) and area-integrated radial Poynting and kinetic flux in the polar outflow at 
r = 190 Rs (right panel). Formally, we define the polar outflow as the region where v r > 
and P < 1. Note the difference in the time range in the panel with the radial fluxes. 

The last panel in Fig. 1 shows the area-integrated radial fluxes of magnetic and kinetic 
energy at r = 190 Rs inside the polar outflow. Formally, we define the polar outflow as the 
region where v r > and f3 < 1. The outflow is Poynting flux-dominated, with the Poynting 
flux exceeding the kinetic energy flux by up to an order of magnitude. 

Figure 2 - Color maps of logarithmic density and toroidal magnetic field overplotted with 
the direction of the poloidal velocity at t — 0.2735 s. The length scale is in units of the BH 
radius (i.e., r' = r/Rs and z' = zj Rs)- 

Figure 3 - Radial profiles of various quantities from our run, time-averaged from 0.2629 
through 0.2818 s. To construct each plot, we averaged the profiles over angle between 
6 = 86° and 94°. The top left panel plots the density (solid line) and temperature (dashed 
line). The top middle panel plots the gas pressure (solid line) and magnetic pressure. The 
top right panel plots the rotational, radial, Keplerian, and Alfven velocities (solid, dashed, 
dot-dashed, and dotted line, respectively), as well as the sound speed (triple-dot dashed 
line). The bottom left panel plots the angular velocity in units of 2c/ R s . The bottom 
middle panel plots the Maxwell stress, a mag , and the Reynolds stress, a gas (solid and dashed 
line, respectively). We calculate the Reynolds stress using eq. (15) in PB03 and show only 
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its amplitude. The bottom right panel plots the radial, latitudinal and toroidal components 
of the magnetic field (dot-dashed, dashed, and solid line, respectively). The length scale is 
in units of the BH radius (i.e., r' = r/Rs). 
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